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Abstract 

In this paper we consider Bayesian estimation for the parameters of 
inverse Gaussian distribution. Our emphasis is on Markov Chain 
Monte Carlo methods. We provide complete implementation of the 
Gibbs sampler algorithm. Assuming an informative prior, Bayes esti- 
mates are computed using the output of the Gibbs sampler and also 
from Lindley's approximation method. Maximum Likelihood and Uni- 
formly Minimum Variance Unbiased estimates are obtained as well. 
We also compute Highest Posterior Density credible intervals, exact 
confidence intervals as well as "percentile" and "percentile-t" bootstrap 
approximations to the exact intervals. A simulation study was con- 
ducted to compare the long-run performance of the various point and 
interval estimation methods considered. One real data illustration 
has been provided which brings out some salient features of sampling- 
based approach to inference. 
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1 Introduction 



The inverse Gaussian (IG) distribution arises as the first p assage time dis- 



tribution in a Brownian motion process with positive drift. iTweedid ( 119571 ) 
first studied its basic characteristics and important statistical properties and 
showed certain analogies between its statistical analysis and that of the nor- 
mal distribution. Chhikara and Folks explored the sampling theory inference 
for the distrib ution in great deta i l over several publications culminating in a 
review paper (IFolks fc Chhikaral . Il978l ). They also proposed it as a lifetime 
model a nd suggested its applicat i on in situations where the initial failure rate 
is high ( jChhikara fc Folka . Il9771 ) . IJohnson et al.l (119941 ) lists many useful ap- 



plications of IG distribution in diverse fields ranging from theoretical physics 
to meteorology, reliability and lifetime data analysis, sequential analysis and 
industrial quality control, business applications and so on. The probability 
density function of the IG distribution is given by: 



f(x\n,X) 



A 



2irx 3 



1/2 



cxp 



X(x — 
2/i 2 x 



x > 0, /i, A > 0, 



which we will denote by lG(fi, A). The mean and the variance of the distribu- 
tion is given by E(X) = fi and Var(X) = fi 3 /X respectively. Hence, neither 
/i can be interpreted as location nor A as scale parameter in the usual sense. 

For a random sample Xi, X 2 , . . . , X n from IG(/i, A), the maximum likeli- 
hood estimators (MLE) of \i and A are given by: 



X and A = n / ^(1/X; - 1/X), 



(2) 



i=i 



where X = ^2 Xj/n, and the uniformly minimum variance unbiased estima- 



i=l 



tors (UMVUE) are given by: 



fi = fi 



X and A = (n -3) /^(1/X - 1/X). 



(3) 



i=l 



Bayesian analysis of the IG distr ibution has received considerable atten- 
tion in the literature. Palmer first considered Bayesian estimation for 
the IG(/i, A) model. He showed that when both parameters are treated as 
unknown, no natural conjugate prior exists and results become very difficult 
to obtain. This led to t he use of other parametriz a tions of IG distribu- 



tion in Bayesian studies. iBanerjee fc Bhattacharvval (119791 ) worked with a 



parametrization involving tp — 1/fM and A and derived some Bayesian results 
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with non-informative reference prior as well as the natural conjugate prior. 
Later works on the (ip, A) parametrization include analysis based on Gibbs 
sampler and a Sampli ng Importanc e Resampling (SIR) scheme, assuming 
non-informative prior ( iShastril . I1994T) and approximate Ba yesian estimation 
using natural conjugate prior ( Ahmad fc Jaheenl . 1995 ). Betro fc Rotondi 
(119911 ) pointed out that in Banerjee and Bhattacharya's approach the pos- 
terior mean of does not exist, that is, the Bayes estimate of the mean 
of the distribution is not available. They chose another parametrization in- 
volving n and (J), where 0~ 1//2 is the coefficient of variation, arguing that it 
renders assessing the priors more natural because of the physical meaning 
attached to the parameters. They considered the marginal prior of to be 
a gamma density and the conditional prior of fi given <p to be an IG density 
and obtained Bayes estimates of the parameters as well as of their inverses. 

Among the works on IG(/i, A) parametrization, iPadgettl (119811 ) consid- 
ered the estimation of reliability function. He noted that the Jeffreys' prior 
tt(h, A) oc (/i 3 A) -1 / 2 does not lead to a proper posterior distribution and using 
the locally uniform non-informative reference prior 7r(/i| A) ex constant and 
7r(A) cx A" 1 found it extremely difficult to compute the e stimat e of rel iability. 



Sinhal (|1986f) solved 



Ismail fc Auda 



As a workaround, he suggested a modified estimator. 
this p roblem by using Lindley's approximation technique. 
(j2006l ) obtained estimates for A and its posterior density when \i is known, 
via Gibbs sampling, assuming Jeffreys' prior for A. 

In this paper, we focus on the IG(/x, A) parametrization and provide full 
Bayesian analysis for the unrestricted model. As a substantial body of work 
is available for the (/z, A) parametrization, especially in the classical frame- 



Palmer ( 


1973); 


Shastri 


( 


1994) 



as possible — a belief also expressed 
compute the Bayes estimates of the parameters under the assumption of in- 
formative priors. If quality prior information is available, then it is preferabl e 
to use informative priors rather than non- informative ones (see iBergerl . Il985[ ) . 

Our emphasis here is on implementation of Markov Chain Monte Carlo 
(MCMC)-based computational methods. We provide complete posterior 
analysis of IG(/i, A) distribution via Gibbs sampler. To put things in perspec- 
tive, we employ different methods of estimation and see how they perform 
compared to each other. Besides computing the Bayes estimates of the pa- 
rameters using the output of a Gibbs sampler, we also compute approximate 
Bayes estimates using Lindley's method. We compare the performance of the 
Bayes estimators with their classical counterparts, namely the MLEs and the 
UMVUEs given by Equation ($2§ and Equation 03]), in terms of their mean 
squared error (MSE) by extensive simulation. Although, the use of frequen- 
tist risk analysis could be suspect from a purely Bayesian point of view, there 
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are convincing argument s in favour of evaluating long run performance of a 



Bayesian procedure (see iBergerl . Il985l . Section 4.8). It is observed that the 



proposed Bayes estimators outperform the classical estimators. 

We also consider different methods of interval estimation and compare 
the performance of the 95% HPD credible intervals — computed from the 
output of the Gibbs sampler - - with the 95% exact confidence intervals 
and two different bootstrap approximations to the exact intervals - the 
"percentile" bootstrap and "percentile-t" bootstrap intervals. Here also, the 
Bayesian credible intervals emerge as better. 

The plan of the rest of the paper is as follows. In Section El we give 
the prior and the posterior distributions. Bayes estimation using Lindley's 
method and Gibbs sampler is outlined in Section[3J SetionHJcontains different 
interval estimation methods discussed in detail. The simulation study is 
described and results are presented in Section [5j In Section El one real life 
dataset is analysed to illustrate our methodology. We summarize our findings 
and conclude the paper in Section [71 



2 Prior and Posterior Distribution 

Let x = (x\,X2, ■ ■ ■ ,x n ) be a random sample of size n from IG(/i, A). The 
likelihood function of the observed sample can be written as: 

/ \ \ n/2 n 

L(x; //, A) = ( — ) J] xT exp [-X(a^ 2 - nfT 1 + 0)] , (4) 
^ ' i=i 

where 

n n 

a = ^x l /2 and ^^x~ l /2 (5) 

1=1 i=l 

and (a, (3) are jointly sufficient for (//, A). 

Under the assumption of independence of /i and A, the joint prior of (/i, A) 
is 7r(/i, A) = 7Ti(/x) 7r 2 (A), where the prior for /i is a gamma density given by 

b a 

7Ti(/i) = — -// a_1 exp [-&//]; a, 6 > 0, (6) 
1 (a) 

and an independent prior for A is given by another gamma density 

vr 2 (A) = ^frX c ~ l exp[-dA]; c, d > 0. (7) 

l (C) 

Then the joint posterior of (/i, A) can be written as 
pO, A| x) = K^- 1 exp[-6/i] A c+ ^ exp [-\(afi~ 2 - n^ 1 + (3 + d)] , (8) 
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where K is the normalizing constant. Therefore, the Bayes estimate of 
h(n,X), any function of fi and A, under squared error loss is the posterior 
expectation hs(fi, A) = J °° h(fi, A) A| x) d/i dA. This integral can not 
be computed in an explicit form. Hence, we discuss two different methods to 
compute the Bayes estimate in the next section. 



3 Bayesian Point Estimation 

In Bayesian Statistics the posterior expectations frequently involve integrals 
which can not be computed explicitly. Various strategies to tackle this prob- 
lem are available, e.g. adaptive quadrature based on numerical analysis, 
analytical approximation methods and sampling-based methods like Monte 
Carlo importance sampling and MCMC methods. Here, we employ one ap- 
proximate Bayes estimation method and one particular MCMC technique to 
compute the Bayes estimates of the parameters. 



3.1 Lindley's Approximation 

Lindley ( 1980l ) suggested an asymptotic approximation to compute the ratio 



of two integrals. Using the prior mentioned in Section EJ the approximate 
Bayes estimates of fi and A under the squared error loss function, given by 
Lindley's method, turn out to be 

fi L = fi+{a + 2)^- h -^r (9) 
nX nX 

and 

o o n A *2cl\ 

X L = X+ 2c -1 , 10 

n n 

where fx and A are the MLEs of /i and A as given in Equation ([2]). We note 
that as n — > oo, fi L — y fi and A^ — > A. The derivations of Equation and 
Equation ffTOl) are given in Appendix A. 

Lindley's method is designed to compute the approximate Bayes esti- 
mates. MCMC methods, e.g. Gibbs sampling, which is discussed next, 
allows us to compute the Bayes estimates as well as to construct HPD cred- 
ible intervals and summarize the marginal posterior distributions effectively 
using samples generated from the posterior. 



3.2 Gibbs Sampling 



Gibbs sampling algorithm was first introduced bv iGeman fc Gemanl fll984l ). 
who applied the sampler on a Gibbs random field. A landmark paper by 
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Gelfand fc Smith! (119901 ) first brought to attention its importance in general 
Bayesian computation. Today Gibbs sampler is the most widely used MCMC 
method and has found application in problems from diverse fields (e.g. for 
ap plication in life testing and reliability as well as details of the algorithm 

see 



Upadhvav et all l20nih . 



From Equation ([S]), the full conditionals for A and fi can be written as 

7r(A| fi, x) oc A c+ ^ _1 exp [-A(a/i~ 2 - n^T 1 + (3 + d)] (11) 
7r(//| A, x) oc fj 1 a ~ 1 exp [— a\fi~ 2 + nXji^ 1 — . (12) 

We see Equation (fTTj) is a gamma distribution and hence A can be easily 
generated using standard available generators. Generating fi from Equa- 
tion (TT21) can be a little tricky as Equation ( TT2|) is not log-concave in \i. We 
provide a procedure in Appendix B, which can be seen as an extension of 
the inverse transform method, which we use to generate samples from Equa- 
tion (|T2|) . This method was tested for several values of the parameters and 
it was found to be quite efficient. 

Now that we have a mechanism to simulate Equation fill I) and Equa- 
tion ( JT2]) . the posterior in Equation flH]) can be easily simulated via a two- 
stage Gibbs sampler. We took single long run of the chain and used ergodic 
averages for convergence monitoring. After discarding the initial transient 
phase, we picked up equally spaced outcomes from the remaining realizations 
of (n, A) and regarded them as samples from the true posterior. The spacing 
was so chosen that it renders the serial correlation negligible. 

Then the Bayes estimates of /i and A under squared error loss function 
are given by 

1 N . 1 N 

fl B = E(fJ,\x) = x^2»j and X B = £(A|x) = Jj^2*3, (13) 

i=i i=i 
where (//,-, Xj) is the j-th MCMC sample, j = 1, 2, . . . , N. 



4 Interval Estimation 

An interval estimate is often more useful than simply a point estimate. Taken 
together, the point estimate and the interval estimate say what is our best 
guess for the parameter value and how far in error that guess might reason- 
ably be. In this section we describe four interval estimation procedures for 

CM). 
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4.1 HPD Intervals 



To construct HPD intervals for u a nd A we follow the Monte Carlo procedure 
proposed bv I Chen fc Shad (ll999l ). Given an MCMC sample (fij,Xj); j = 
1, 2, . . . , N, we calculate the HPD interval for \i as follows: 
Step 1. Sort {fij, j = 1,2, ... , N} to obtain the ordered values 

/^(l) < H(2) < ■ ■ < H(N)- 

Step 2. Compute the 100(1 — a)% credible intervals 

Ri(N) = (//(i), /i( i+[ (i_ a )iv])), for i = l,2,...,N-[(l-a)N], 
where [(1 — a)N] is the integer part of (1 — a)N. 

Step 3. Ri*(N), the 100(1 -a)% HPD interval, is the one with the smallest 
interval width among all _Rj(iV)'s. 

The same procedure can be applied to calculate the HPD interval for A. 



4.2 Exact Confidence Intervals 



Confidence intervals for A parallel those for the variance of the normal dis- 
tribution. A 100(1 — a)% two sided confidence interval for A follows from 



the fact that n\V ~ xL-i)> where V = Hi l / X i ~ l / x )/ n = V A (ITweedid . 



19571 ) . and is given by 



i=l 



^n~l,a/2 ^n-l,l-a/2 



nV 



nV 



(14) 



where Xn-i a * s ^e (100 ■ a)th percentile of x 2 distribution with {n — 1) 
degrees of freedom. 

For un known A, a 100(1 — a)% two sided confidence interval for ji was 
derived by IChhikara fc Folks! (119761 ) and is given by 



X[l + JXV/(n - 1) ti_ Q/2 ]-\ X[l- JXV/(n - 1) ti- a/2 



i-i 



(15) 



if (l - y/XV/(n - 1) t!_ Q/2 ) > 0, and (x[l + ^XV/(n - 1) t^ a/2 ]-\ oo 

otherwise, where tx-a.li is the 100(1 — a/2)th percentile of Student's t distri- 
bution with (n — 1) degrees of freedom. 
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4.3 Bootstrap Confidence Intervals 

In this subsection we discuss two bootstrap methods for constructing confi- 
dence intervals. We use the "percentile" or bootstrap-p method, which is the 
most widely used one and the "percentil e-t"or bootstrap-t metho d, which is 
the most efficient one theoretically (see lEfronl . Il982l ; lHalll . Il988l ). Our con- 
structions are based on the parametric bootstrap. We describe the method- 
ology for fi. Intervals for A can be constructed in a similar way. 

The bootstrap-p confidence interval for \i can be obtained as follows: 
Step 1. Calculate ft and A from x. 

Step 2. Generate B bootstrap samples x* of size n from IG(jl, X) and 
calculate fi* and A* from x*, i — 1, 2, . . . , B. 

Step 3. The 100(1 — a)% bootstrap-p confidence interval for is given by 



~*(ct/2) ~*(l-a/2) 



where is the (B ■ a)th value in the ordered list of the B replications of 
/'•• 

The bootstrap-t confidence interval for /i can be obtained as follows: 
Step 1. Calculate fi and A from x. 

Step 2. Generate a bootstrap sample x* of size n from IG(£l, A) and calculate 
/}* and A* from that. 

Step 3. Generate a bootstrap sample x** of size n from IG(fi*,X*) and 
calculate ft* and A** from that. 
Step 4. Repeat Step 3 B 2 times. 

/ ~B~2 B 2 

Step 5. Calculate se(/2*) = J^- E(£f ~ fi**) 2 , where fi** = J2 fi7/B 2 - 

V 2 i=i j=i 

Step 6. Calculate T* = ^ ~ ^ . 

Step 7. Repeat Steps 2-6 Bi times. 



/ B 2 B 2 

Step 8. Calculate se{p) = W^r Yifij ~ fi*) 2 ' where fi* = Yl fij/ B 2- 

V 2 3=1 3=1 

Step 9. The 100(1 — a)% bootstrap-t confidence interval for \x is given by 
{fi - T* B f~ a/2) ■ se(fi), fi - T B [ a/2) ■ se(fi)) , 

where T B a ^ is the (Bi ■ a)th value in the ordered list of the B\ replications 
ofT*. 
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5 Simulation Study 



To compare the performance of the various estimation methods described 
above, a simulation study was performed. We took samples of size n = 
15,20,30 and 50 from IG(3,4). We assumed gamma(6,2) as prior for \i 
and gamma(5, 1.25) as prior for A, i.e. we chose a = 6,b = 2,c = 5 and 
d = 1.25. Under these priors, the Bayes estimates under squared error loss 
for fi and A were computed for each sample using both Lindley's method and 
the output of the Gibbs sampler. MCMC samples of size 1000 were taken 
for the computation. We also computed the MLEs and the UMVUEs. The 
average value of the estimates and their MSE, based on 1000 replications, 
are reported in Table [TJ 

The 95% exact, bootstrap-p and bootstrap-t confidence intervals and 
HPD intervals for the parameters were also computed for each sample. For 
each interval given by (I, u), we evaluated a shape factor measuring the skew- 
ness of the interval, which is (u — MLE)/(MLE — I) for the confidence intervals 
and (u — posterior mode)/(posterior mode — I) for the HPD intervals. We 
also calculated two factors as MissLeft and MissRight for each interval, where 
MissLeft implies that the true value of the parameter falls below I and Mis- 
sRight means that it lies above u. Table |2] gives the average intervals along 
with their average shape and coverage, including the MissLeft and MissRight 
probabilities, based on 1000 replications. 



Table 1: Average and MSE (within parenthesis) of different estimators for 
different sample sizes. The first entry in each cell corresponds to fi and the 
second to A. 



n 


MLE 


UMVUE 


Bayes(Lindley) 


Bayes(Gibbs) 


15 


3.0028(0.4396) 
4.9167(5.4989) 


3.0028(0.4396) 
3.9333(2.9860) 


3.2103(0.3444) 
3.0613(12.3134) 


3.1426(0.2764) 
4.1909(0.8051) 


20 


3.0087(0.3212) 
4.6511(3.2504) 


3.0087(0.3212) 
3.9535(2.0442) 


3.1831(0.2799) 
3.6867(0.7075) 


3.1357(0.2343) 
4.1771(0.7859) 


30 


3.0076(0.2220) 
4.4509(1.6967) 


3.0076(0.2220) 
4.0059(1.2096) 


3.1327(0.2072) 
4.0109(0.3317) 


3.1106(0.1868) 
4.1758(0.6725) 


50 


2.9924(0.1363) 
4.2756(0.8509) 


2.9924(0.1363) 
4.0191(0.6851) 


3.0734(0.1312) 
4.0924(0.4033) 


3.0663(0.1260) 
4.1363(0.4879) 



Table [T] brings certain points to our attention. First of all, it is observed 
that as expected, the MSE of each estimator decreases as the sample size 
increases. Bayes estimators obtained from both Lindley's method and Gibbs 
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Table 2: Average 95% intervals and their shape and coverage probabilities 
for different sample sizes. The first entry in each cell corresponds to // and 
the second to A. 



n 




Interval 


Shape 


Coverage 


MissLeft 


MissRight 




Exact 


(2.0359, 6.0797) 
(1.8450, 8.5612) 


2.9480 
1.1865 


0.956 
0.960 


0.018 
0.021 


0.026 
0.019 


15 


Boot-p 


(1.9349, 4.4723) 
(2.8191, 13.1053) 


1.3519 
3.9045 


0.911 
0.862 


0.006 
0.138 


0.083 
0.000 




Boot-t 


(2.0587, 5.4791) 
(1.4864, 8.8460) 


2.5083 
1.1465 


0.944 
0.952 


0.023 
0.025 


0.033 
0.023 




HPD 


(2.0346, 4.4732) 
(2.0341, 6.5653) 


1.9159 
1.5038 


0.973 
0.989 


0.008 
0.000 


0.019 
0.011 




Exact 


(2.1441, 5.1679) 
(2.0713, 7.6400) 


2.3929 
1.1585 


0.956 
0.955 


0.020 
0.021 


0.024 
0.024 


20 


Boot-p 


(2.0507, 4.2688) 
(2.8284, 10.4286) 


1.3004 
3.1769 


0.927 
0.896 


0.008 
0.104 


0.065 
0.000 




Boot-t 


(2.1707, 4.9529) 
(1.8732, 7.7919) 


2.2531 
1.1303 


0.950 
0.953 


0.023 
0.026 


0.027 
0.021 




HPD 


(2.1334, 4.3196) 
(2.1918, 6.3284) 


1.8339 
1.4441 


0.974 
0.981 


0.008 
0.002 


0.018 
0.017 




Exact 


(2.2764, 4.4662) 
(2.3808, 6.7836) 


1.9528 
1.1268 


0.953 
0.951 


0.020 
0.032 


0.027 
0.017 


30 


Boot-p 


(2.1987, 4.0202) 
(2.9200, 8.3060) 


1.2437 
2.5194 


0.930 
0.905 


0.010 
0.095 


0.060 
0.000 




Boot-t 


(2.2905, 4.4334) 
(2.2752, 6.8739) 


1.9542 
1.1154 


0.943 
0.953 


0.025 
0.026 


0.032 
0.021 




HPD 


(2.2586, 4.0954) 
(2.4484, 6.0391) 


1.7282 
1.3764 


0.961 
0.981 


0.013 
0.004 


0.026 
0.015 




Exact 


(2.4047, 3.9693) 
(2.6983, 6.0049) 


1.6465 
1.0964 


0.951 
0.947 


0.019 
0.030 


0.030 
0.023 


50 


Boot-p 


(2.3468, 3.7587) 
(3.0425, 6.7622) 


1.1836 
2.0179 


0.932 
0.921 


0.009 
0.077 


0.059 
0.002 




Boot-t 


(2.4133, 3.9833) 
(2.6452, 6.0585) 


1.6970 
1.0947 


0.939 
0.954 


0.022 
0.027 


0.039 
0.019 




HPD 


(2.3877, 3.8254) 
(2.7222, 5.6353) 


1.5774 
1.2881 


0.956 
0.968 


0.012 
0.012 


0.032 
0.020 
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sampler are observed to perform much better than the classical estimators, 
but the discrepancy in their relative performance tends to get smaller and 
smaller with the increase in sample size. We observe that except for estima- 
tion of A in small samples, Bayes estimators obtained from Lindley's method 
perform quite well. This can be explained by the rather erratic behaviour 
of the MLE of A in small samples. As the performance of the MLE of A 
gets stabilised with larger sample size, Lindley's method quickly catches up 
and performs much better for moderate to large sized samples. Overall, the 
Bayes estimators computed from our proposed Gibbs sampler are found to 
be performing very well. 

From Table [2] we see that as par our expectation, length of all the inter- 
vals decrease with increase in sample size. We also observe that the cover- 
age probability of each interval remains quite close to the nominal value of 
95%. As the exact confidence interval, and hence the bootstrap intervals, 
are equi-tailed by construction, the ideal values of MissLeft and MissRight 
probabilities are 0.025 each. The bootstrap-t approximates the shape of the 
exact intervals very well and as a result gets both the endpoints and balance 
at both ends right. The bootstrap-p gets the tail area probabilities wrong for 
both MLEs and hence, as approximation of exact intervals, it fares poorly in 
comparison to bootstrap-t. It is seen that for all the intervals including the 
HPD, shape — > 1 with increase in n, i.e. the intervals become more and more 
symmetric. This results from the asymptotic normality of the underlying 
distributions of the MLEs and the posteriors of fi and A. For symmetrical 
posteriors, the HPD interval becomes equi-tailed. Indeed, we observe that as 
sample size increases, the MissLeft and MissRight probabilities for the HPD 
get stabilised toward 0.025. Overall, we see that the HPD intervals work 
very well in terms of both length and coverage. 



6 Numerical Illustration 



In this section, we consider a real data set initially analysed by lChhikara fc Folks 
fll977t ). The maintenance data set given in Table [3] represents active repair 
times (in hours) for an airborne co mmunication transceiver . Among the 
previous works done on this data set, IChhikara &: Folks! (119771 ) obtained un- 
biased and maximum likelihood estimation of IG parameters and showed 
that the I G model is a good fit to the data using Kolmogorov-Smirnov test. 
Sinhal (119861 ) computed the Bayes estima t es usi ng Lindley's method under a 
non- informative prior. iBetro fc Rotondil ( )199ll ) also analysed this data set 
and calculated the estimates for the parameters. 

For our analysis, in absence of prior information, we assumed the non- 
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Table 3: Repair times (in hours) of 46 transceivers 



0.2 0.3 

0.7 0.8 

1.5 1.5 

3.3 3.3 

7.5 8.8 



0.5 0.5 

0.8 1.0 

1.5 2.0 

4.0 4.0 

9.0 10.3 



0.5 0.5 

1.0 1.0 

2.0 2.2 

4.5 4.7 

22.0 24.5 



0.6 0.6 

1.0 1.1 

2.5 2.7 

5.0 5.4 



0.7 0.7 

1.3 1.5 
3.0 3.0 

5.4 7.0 



in formative prior 7r(/x| A) oc constan t and 7r(A) oc A^ 1 , the same prior used 
bv IPadgettl (I198ll ) and ISinhal (119861 ). The joint posterior of (/i, A) under this 
prior can be written as 



p(/i, A| x) oc A 2 1 exp [— X(ajd~ 



n/i 



1 +fi)]. 



(16) 



This posterior can be obtained by putting a=l,5 = 0,c = and d = in 
Equation (jHJ). So, for our analysis we could use the same Gibbs sampler as 
described in Section I3~2"} without any further modification, for this particular 
choice of a, b, c and d. We used MLE as the initial value for the chain and took 
MCMC samples of size 1000 to compute the Bayes estimates. Some sample 
based posterior characteristics are presented in Table HI Figured] depicts the 
marginal posterior density estimates of \i and A using the Gaussian kernel, 
based on the MCMC samples. The interval estimates obtained, are presented 
in Table EJ 



Table 4: Sample based posterior characteristics for the repair times data 





Mean 


Variance 


Mode 


1st Quartile 


Median 


3rd Quartile 


Minimum 


Maximum 


/' 


4.3999 


1.8245 


3.7475 


3.4750 


4.0229 


4.9173 


2.2434 


9.9877 


A 


1.6129 


0.1164 


1.6065 


1.3791 


1.5953 


1.8092 


0.7234 


2.9163 



Using the vague prior mentioned above, ISinhal (119861 ) obtained the ex- 
pressions for the Bayes estimates of /i and A, using Lindley's method, as: 



\x L = ^ H and \ L 

nX 



n — 1 



n 



A. 



(17) 



and Equation fTTTJ]) for the 



which, again, can be derived from Equation 
same choice of a, b, c and d as above. 

The MLEs, UMVUEs and Bayes e stima tes from Lindley's method for 
the parameters, as obtained by ISinhal (119861 ). are £i = 3.6065, A = 1.6589; 
fj, = 3.6065, A = 1.5507; fi L = 4.1178 and X L = 1.6228 respectively. The 
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Table 5: 95% intervals and their shape for the repair times data 



Exact Boot-p Boot-t HPD 

Interval (2.4998, 6.4715) (2.2868, 5.4372) (2.4251, 6.5312) (2.4134, 7.2643) 

Shape 2.5888 1.3872 2.4756 2.6361 

Interval (1.0229, 2.3588) (1.1950, 2.7510) (0.8674, 2.4426) (0.9268, 2.2686) 

Shape 1.1007 2.3547 0.9903 0.9742 



Bayes estimates obtained from our proposed Gibbs sampler turn out to 
be ur = 4.3999 and A r = 1.6129. The Bayes estimate of \i provided by 



Betro fc Rotondil (Il99ll ) is 3.653. So, the estimates obtained are in reason- 
able agreement with previously proposed estimates. 



>^0.3 




1.0- 



-go.5- 



0.0- 




lambda 



Figure 1: Kernel density estimates for marginal posterior densities of \i and 
A for the repair times data 

Figure [1] and Table H] reveals the advantage of the MCMC-based or 
sampling-based approach towards inference. Full posterior analysis based 
on MCMC samples, gives an overall view of the underlying posteriors and 
suggests the use of appropriate estimators for inference purpose. For exam- 
ple, for the repair times data, unlike in the case of A, a little difference can be 
seen between possible estimates of \i. Now, Figure [T] shows that the marginal 
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5 10152025 5 10152025 5 10152025 



Figure 2: The empirical and fitted CDFs for the repair times data using 
different estimates: left: MLE for \i and A centre: posterior mean for \i and 
A right: posterior mode for fi and posterior mean for A 

posterior of \i is much more skewed than that of A. This suggests that for 
/i, HPD intervals are more preferable than equi-tailed intervals and modal 
value might be a more reasonable point estimate. Indeed, from Figure |5] we 
can see that using posterior mode instead of posterior mean for /i, we can 
obtain a better fit to the data. 

7 Conclusion 

In this paper, we have considered Bayesian estimation for the parameters 
of the IG(/i, A) distribution. We assumed independent gamma priors for 
both the parameters in this study. Complete implementation of the Gibbs 
sampling algorithm has been provided. Simply by choosing suitable hyper- 
parameters, and without any further modification, the proposed Gibbs sam- 
pler is shown to work as smoothly for a particular non-informative prior, 
which has been used previously for studying the IG distribution. We have 
used squared error loss function in this paper, but the proposed methodol- 
ogy can be routinely applied for other asymmetric loss functions as well. We 
have shown how techniques based on MCMC can easily deal with the is- 
sue of awkward posterior, perform well compared to other available methods 
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of estimation and on top of that, provide a conceptually simple, routinely 
implementable and unified framework for complete analysis of the IG distri- 
bution. But perhaps the greatest advantage of this sampling-based approach 
is the unique insight it offers us into the posteriors. This enables us to make 
better inference by helping us choose appropriate inference procedures, in 
particular when the posterior turns out to be non-normal and skewed. 



Appendix A 

For a model involving two parameters (9\, 6 2 ), Lindley's approximation tech- 
nique ( Lindleyl . 1980| ) gives an expression for the posterior expectation of u, 



a function of (^i,^), as: 



2 2 



E(u\x) = u + - ^ ^2 u » a ^ + U iPi + ^ L ni a n u i 

i=l j=l j=l 

+ ^Lu 2 (2a 12 Ui + a u U 2 ) + -^122(0-22^1 + 2cr 12 U 2 ) 

+ ^L 222 a 22 U 2 , (18) 

where all the terms are evaluated at the MLE (#1, 8 2 ) and Uk = Yl u i a ki, t ne 

i 

partial derivatives pj, Uj, Uij, Lyfc are defined as, for example, = , 
where each suffix denotes differentiation once with respect to the variable 
having that suffix, are the elements of the inverse of the matrix having 
elements {— L^}, L is the log-likelihood and p is the logarithm of the joint 
prior. 

Now, considering (9i,9 2 ) = (//, A), the definitions of likelihood function 
and priors given in Section [2] lead to: 

n n 
- — ,Li 22 = 0,L 222 = j-, 
^ A 3 

n 2X2 
1 = U, a 22 = , 

n 









Lin 




, -^112 




£< 












= 4 




cr n = 




0"12 = 




nA 






a — 1 




Pi = 




-b 









A 

Considering u = p and u = A in turn and using the above values, from Equa- 
tion (fl8l) we get the expressions of Bayes estimators given by Equation Q 
and Equation (TlOT) respectively. 
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Appendix B 



Let X be a random variable having probability density function /(■) and 
cumulative distribution function F(-). Then F(x) ~ {7(0, 1), and this result 
is the foundation of the inverse transform method of random variate genera- 
tion. We generate a u from {7(0, 1), and find out c such that F(c) = u. Then 
c can be taken SIS db random sample from /(•). 

Now, let f(x) = kg(x), where k is a constant, i.e. g(-) is the non- 
normalized density. Then, defining F*(x) = f* g{t)dt = F(x)/k, we have 
F*(x) ~ {7(0, 1/k). So, generating u from {7(0, 1/k) and solving for c in 
F*(c) = u, also gives us a random sample c from /(•). 

We use this idea to generate /i from Equation (IT2|) . First, we find a c 
such that J C 7r(/x|A, x)d/x ~ Jq 00 7r(//|A, x)d/x, i.e. integrating 7r(/x|A,x) upto 
c, yields almost — for example, 99.9999% of — the whole area under the 
curve. Let J" ° 7r(/x|A, x)d/i = T. We generate a u from {7(0, T) and solve 

Jq it(ji\X, x)d/x = m for /i*. Then /i* is a random sample from 7r(/x|A,x). 
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